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Summary 

We proposed a novel characterization of errors for numerical weather predictions. A gen- 
eral distortion representation allows for the displacement and amplification or bias correction 
of forecast anomalies. 

Characterizing and decomposing forecast error in this way has several important appli- 
cations, including the model assessment application, and the objective analysis application. 
In this project we have focused on the assessment application, restricted to a realistic but 
univariate 2-dimensional situation. Specifically we study the forecast errors of the sea level 
pressure (SLP), the 500 hPa geopotential height, and the 315 K potential vorticity fields for 
forecasts of the short and medium range. The forecasts are generated by the GEOS (Goddard 
Earth Observing System) data assimilation system with and without ERS-1 scatterometer 
data. 

A great deal of novel work has been accomplished under the current contract. In broad 
terms we have developed and tested an efficient algorithm for determining distortions. The 
algorithm and constraints are now ready for application to larger data sets, to be used to 
determine the statistics of the distortion as outlined above, and to be applied in data analysis 
by using GOES water vapor imagery to correct short term forecast errors. 
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We proposed a novel characterization of errors for numerical weather predictions (Hoffman 
et al. 1995 [7]). In its simplest form we decompose the error into a part attributable to 
phase errors and a remainder. The phase error is represented in the same fashion as a 
velocity field and is required to vary slowly and smoothly with position. A general distortion 
representation allows for the displacement and amplification or bias correction of forecast 
anomalies (Hoffman and Grassotti 1996 [6]). 

Characterizing and decomposing forecast error in this way has several important ap- 
plications. For the model assessment application, our approach results in new objective 
measures of forecast skill which are more in line with subjective measures of forecast skill 
and w'hich are useful in validating models and diagnosing their shortcomings. With regard 
to the objective analysis application, meteorological analysis schemes balance forecast error 
and observational error to obtain an optimal analysis. Presently, representations of the error 
covariance matrix used to measure the forecast error are severely limited. For the objective 
analysis application our approach will improve analyses by providing a more realistic mea- 
sure of the forecast error. We expect, a priori , that our approach should greatly improve the 
utility of remotely sensed data which have relatively high horizontal resolution, but which 
are indirectly related to the conventional atmospheric variables (Hoffman and Grassotti 1996 
[6]). A related application is to combine two sources of data, e.g. SSM/I and ground-based 
radar rain rates (Grassotti et al. 1998 [4]). 

In this project we have focused on the assessment application, restricted to a realistic 
but univariate 2 dimensional situation. Specifically we study the forecast errors of the sea 
level pressure (SLP), the 500 hPa geopotential height, and the 315 K potential vorticity 
fields for forecasts of the short and medium range. The forecasts are generated by the 
GEOS (Goddard Earth Observing System) data assimilation system with and without ERS- 
1 scatterometer data. Our studies are a first step towards (1) a testbed for the use of 
the distortion representation of forecast errors, (2) a means of validating the GEOS data 
assimilation system and (3) a description of the impact of the ERS 1 scatterometer data. 

This final report for NAS5-32953 is also offered in lieu of the annual report for the third 
and final year of the project. We report on the entire research project in this report. The 
particular areas of research studied during this past year are: 

1. Extending the analysis to potential vorticity at 315 K. Examples are given of the 
March 1993 super storm. 

2. Allowing for time continuity by treating the increment in C as the control variable to 
be constrained. (The C are the spectral representation of the distortion.) 

3. Developing an approach to determine the statistics of the C . 

4. Experimenting with different stopping criterion for the stepwise truncation part of our 
approach to determine the statistics of the C, in terms of isostropv and homogeniety. 



2 


AER, Inc. P599, Final Report 


5. Aquiring ECMWF Lorenz data sets of 500 hPa height. These data sets conviently 
provide large samples for determining the statistics of the C. 

6. Developing an approach to use 6.7 pm water vapor imagery. To quality control cloud 
contamination we will use a window channel. 

7. Aquiring OPTRAN and coefficients for GOES. GOES data are acquired routinely by 
AER for the northern hemisphere. 

2 Data 

Forecasts and verifying analyses made with the GEOS data assimilation and forecast system 
(Schubert et al. 1993 [20]) are used here. The particular experiments studied are described 
by Atlas et al. (1995 [1]) in a study of the impact of ERS-1 scatterometer data on numerical 
weather prediction. The period of study is March, 1993. The forecast model and data 
assimilation system used in these experiments are identical to the GEOS-1 system described 
by Schubert et al., except for some minor bug fixes and the modifications necessary to 
utilize surface wind vectors. Thus the control forecasts in the impact study are standard 
GEOS forecasts. In addition to the CONTROL experiment, several using different types of 
scatterometer wind information are available. Our initial prototyping and sensitivity studies 
use only the 2 x 2.5° CONTROL forecast for the period 6-11 March 1993. In addition 
we have made some comparisons to the corresponding PGLA and VARGLA forecasts. In 
all cases the CONTROL GEOS data assimilation is used as verification. The PGLA and 
VARGLA forecasts use the same setup as the CONTROL forecast, but both add ERS 1 
scatterometer data to the CONTROL data sets in determining the initial conditions for the 
forecast. In PGLA the scatterometer data is processed using the directional filtering method 
of Offiler (1994 [12]), while in VARGLA, the variational analysis method of Hoffman (1984 
[5]) is applied to the ERS-1 o° measurements. 

3 Methodology 

Since we require that these distortion fields vary smoothly, a spectral representation is ap- 
propriate. Determining the distortion which provides the best match is then equivalent to 
minimizing the misfit between the first field and a distortion of the second, w r ith respect to the 
spectral coefficients of the distortion. In the present project we use a global or hemispheric 
domain, and spherical harmonics as basis functions. 

In brief, the distortion is determined by minimizing the objective function J, by varying 
the displacement arid bias correction fields, where 

3 — 3 j. -j- 3d -j- 3d. 

The residual cost function, 3 r , measures the misfit of the distorted forecast to the verifying 
analysis. Minimizing 3 r improves the agreement between the (distorted) forecast and the 
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analysis. The two additional penalty terms in the objective function, J d and J a , ensure 
that the final distortion produced by the minimization is relatively smooth and not too 
large. (The terms cost function, objective function and penalty function are used more or 
less interchangeably in the literature. Here, the objective function is the quantity to be 
minimized, a cost function measures lack of fit to data and a penalty function measures lack 
of fit to a constraint.) The smoothness penalty function, Jd measures the roughness of the 
x- and ^-displacements and of the bias correction, ensuring that the distortion is large scale. 
The barrier penalty function, J a , measures the magnitude of the distortion components in a 
way so that small distortions are not penalized, but large distortions are penalized heavily. 
This has the effect of setting up a barrier to the size of the distortions which are determined. 
These last two terms are evaluated using the spectral coefficients of the distortion. 

The three terms making up J are described in the following sections. However, in our 
work so far, the spectral truncations used are so severe that J d and J a are not used in 
obtaining the results presented here. In our studies at the beginning of the project we found 
that using the barrier and smoothness penalty functions results in distortions which are 
smaller in magnitude, but larger in scale, and residual errors which are larger in scale and 
magnitude. The spectra of the original and residual forecast error, both with and without the 
penalty functions show that a great deal of the forecast error on the scales of the distortion 
is explained by the distortion and that the penalty constraints have a strong effect limiting 
the smallest scales in the distortion. 

Note that the limits used to define J a are found to be very useful to precondition the 
minimization, even in cases where J a is not used in the functional. Simulation experiments 
demonstrated that if the control vector is scaled by its limiting values estimate, the true 
solution is quickly recovered. If the scaling is derived from the smoothing function instead, 
the minimization quickly fails with false convergence. For the case of uniform scaling of the 
control vector, the minimization is only partially successful: the objective function is reduced 
only slowly, and after 100 iterations, only half of the original forecast error is explained. 

3.1 Residual cost function, J r 

The residual cost function J r measures the misfit between the distorted forecast and the 
verifying analysis. We denote the forecast by F, the distorted forecast by P, and the verifying 
analysis, or what is considered truth, by T. The cost function is 

UP - Tfda 
T ~ I'do ’ 

where the integral is the surface integral over the global domain. The distorted forecast P 
is obtained from the unmodified forecast F by adding a location-dependent bias correction 
B(X, 6) to the values displaced by the displacement vector field D(A, 9) = ( D U ,D V ), where 
D is expressed here in terms of its zonal and meridional components, in analogy to a wind 
field. Thus, we may write 


P(\,0) = F{\',6 , ) + B(\,0) 
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where the location (A', O') is found by following the displacement vector D( A, 0) back from 
its endpoint (A, 0). 

We represent the scalar field B by a truncated series of spherical harmonics, and the vector 
field D in terms of the spectral coefficients of the corresponding vorticity (C) and divergence 
(5) fields. A degree of smoothness can thus easily be imposed by the truncation of the series, 
and further constraints can separately be imposed on the divergent and rotational parts of 
the displacement field. The control vector C for the optimization problem is thus composed 
of the spectral coefficients for B, (, and 5: 

C = (B, (,S) T . 

Both the forecast F and the verifying analysis T are available on regular latitude- 
longitude grids. For evaluation of the integral, it is convenient to first interpolate T to 
a Gaussian latitude-longitude grid, in which case the formula for J r takes the form 

j 3 i 

where indices i, j denote the grid point location in longitude and latitude, Nj is the number 
of longitude points for latitude j (this number will depend on j only for reduced Gaussian 
grids), and wj is the Gaussian weight for latitude j. These weights are normalized such that 
their sum over all latitudes is unity. 

The first step in the evaluation of the requires the spectral transformation from C 
to B i:j and (D tt ,A,) y . The next step is the evaluation of F(X',0'). Following Ritchie (1987 
[18]), we define latitude-longitude points in terms of 3-dimensional cartesian vectors centered 
on the unit sphere. The origin point (A', O'), corresponding to the 3-dimensional cartesian 
vector r, is then found in the plane of the endpoint location vector g (corresponding to 
gridpoint (A;, Oj)), and the displacement vector d (corresponding to (D u , D v )ij): 

r = ag + /3d, 

where the coefficients a and j3 are chosen to satisfy the constraint that r must lie on the 
surface of the sphere, and that the length of the displacement vector d is equal to the 
great circle distance between g and r. Finally, the value F(X',0') is obtained by bilinear 
interpolation in longitude and latitude from the surrounding grid point values. 

3.2 Smoothness penalty function, Jd 

The smoothness penalty function, J d is given by a simple quadratic form in terms of the 
spectral coefficients of the distortion, 

J d = T>jWj(Cj - Cj) 2 , (1) 

where j ranges over the ordering of the spectral wavenumber vectors, k, and over the com- 
ponents of the distortion— J3, C, S. The C are the background or first guess of the distortion. 
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Thus it. is the increments of C which are constrained. A reasonable first guess is the estimate 
of C obtained for the same forecast, but valid 6 or 12 h earlier. This is a method to enforce 
consistency in the evolution of the distortion. Consistency is most lacking if a feature in 
the forecast may be explained in terms of more than one feature in the verification. These 
situations are more likely to occur as forecast length increases. Different specifications of Wj 
are possible. For example, consider the part of Jd due to the bias correction. In continuous 
form this is given by, 

J6,b = 4 / /(V-[B - B]) 2 - 

Here v is an adjustable parameter normally taken to be 1 and 0 b is the scale for B. Larger 
values of u result in greater smoothing by emphasizing the contributions of higher wavenum- 
bers to Jd- Using the spectral representation of B , the orthonormality of the spherical 
harmonics \F„ m , and the eigenstructure V 2 T„ m = —n(n + l)^ n m , we find that, 

Wk,B = (n k (n k + \)) 2v /o 2 B . 


3.3 Barrier penalty function, J a 

The functional J a serves to limit the amplitude of the distortion. For efficiency the limits 
are set on the spectral coefficients. These limits are chosen in such a way that the grid 
point, (or physical space) values of bias and displacement at all locations will be limited 
by specified values. In addition the spectral coefficient limits provide a good scaling (or 
conditioning) for the minimization. Test runs using synthetic data indicate that convergence 
of the minimization is sensitive to the scaling of the control vector. 

The form of J a is chosen to be, 

Ja = £>(!<?, - Ql/S,) 2 ", 

where Sj are the spectral limiting values for the increment of Cj. The adjustable parameter 
//, nominally 10, controls the steepness of the barrier in spectral space (Fig. 1). 

There is no unique way of setting the spectral limits. We choose limits which correspond 
to an equipartitioning, among the spectral modes, of the contributions to the physical space 
bias correction or displacement component. Here mode means each pair (m, n). The rea- 
soning for this is that no matter what the signs of the spectral coefficients, the modes will 
tend to add up somewhere in the physical domain. On the other hand, the contributions 
within a particular mode, for example due to the sine and cosine components, are always 
out of phase and therefore add up in an rms sense. The limits on components are chosen to 
correspond to a further equipartitioning. 

3.4 Implementation details 

The algorithm is implemented in Splus and Fortran. The spectral transform and computation 
of Gaussian latitudes and weights use a set of general purpose Fortran library functions. All 
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Figure 1: The component of the barrier function for a single term (y = x where x — Cj/ Sj), 
for p equal to 1, 2, 10, 20. 


computations are performed in double precision. To minimize J we use the built-in Splus 
function nlminb, which implements the algorithms of Dennis et al (1981 [3]). The nlminb 
algorithm uses function and gradient values. Second derivatives of the cost function are 
estimated by finite differences, using repeated evaluations of the gradient and cost function. 

At first we used the version of the minimization algorithm which requires function values 
only. However, the finite difference approach is computationally inefficient, and we have 
recently developed the adjoint of the calculation of J . To develop the adjoint we use tools 
previously developed for this purpose (Hoffman et al. 1992 [8]). In addition, in the present 
case the spectral transforms are nearly self-adjoint (Hoffman and Nehrkorn 1989 [9]), so 
the amount of actual adjoint code for the transforms is limited. The adjoint calculates 
the gradient of the cost function very efficiently. This technique gains more than an order 
of magnitude decrease in computational time and provides a more accurate gradient. The 
efficiency of the adjoint code allows us to use more than enough iterates. All results presented 
here use 100 iterations. 

The minimization starts from the reasonable initial estimate of zero distortion. This 
is a point of maximum non-differentiability. To eliminate this problem, we interpolate the 
original analysis to a new grid using Gaussian latitudes and longitudes offset by half a grid 
length. Now a zero displacement corresponds to locations interior to the grid of the analysis, 
where the interpolation of the analysis is differentiable. 
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Results obtained from our method during this past year extend the results reported previ- 
ously to additional variables, forecasts, and truncation parameters of the distortion. Sen- 
sitivity tests involving different choices of truncation parameters are described in Section 
5.2. We report here on the results of distortion experiments for isentropic potential vorticity 
(IPV), concentrating on the 315 K isentropic surface (q315), and compare those with results 
for the 500 hPa geopotential height (h500) and sea level pressure (sip). We have computed 
distortions for all three variables, for the CONTROL forecast started at 00 LTC 11 March 
1993, for forecasts from 12 to 120 hours, every 12 hours. This forecast period includes the 
so-called ” Superstorm 93”, an exceptionally intense cyclone on the East Coast of the United 
States. Truncation sensitivity tests were also performed for the CONTROL forecast started 
at 00 UTC 6 March 1993 (see Section 5.2). 

A summary of the minimization results is shown in Table 1, which shows the statistics 
of the minimization. 


Table 1: Summary statistics for 315 K IPV (q315), 500 hPa height (h500), and sea level 
pressure (sip) for the CONTROL forecast every 12 h. 


Forecast 

Hour 

sip 

Rms Error 

Var. 

Expl. 

h500 
Rms Error 

Var. 

Expl. 

q315 

Rms Error 
— 

Var. 

Expl. 

Initial 

Final 

Initial 

Final 

Initial 

Final 

012 

1.495 

0.761 

74.1 

20.323 

8.309 

83.3 

0.296 

0.227 

41.1 

024 

2.383 

1.049 

80.6 

26.775 

10.298 

85.2 

0.375 

0.271 

47.8 

036 

3.549 

1.521 

81.6 

33.514 

11.247 

88.7 

0.465 

0.331 

49.3 

048 

4.903 

1.993 

83.5 

42.399 

13.494 

89.9 

0.601 

0.362 

63.7 

060 

6.005 

2.084 

88.0 

51.822 

15.881 

90.6 

0.664 

0.378 

67.6 

072 

6.627 

2.247 

88.5 

62.590 

14.795 

94.4 

0.770 

0.414 

71.1 

084 

7.473 

2.298 

90.5 

73.342 

17.407 

94.4 

0.862 

0.417 

76.5 

096 

8.577 

2.296 

92.8 

86.126 

21.365 

93.8 

0.933 

0.445 

77.2 

108 

9.748 

2.436 

93.8 

103.189 

24.016 

94.6 

0.965 

0.483 

75.0 

120 

11.176 

2.512 

94.9 

118.070 

27.080 

94.7 

1.021 

0.504 

75.6 


In Table 1 the rms error is the rms residual error for the hemisphere (in m, hPa , and 
PV units, respectively), calculated as the square root of J T . Initial values are for a zero 
distortion, final values for the solution. The last column is the fraction of the initial error 
variance (J r ) explained by the distortion (in %). 

In the reported cases 100 iterations were used by the minimization, and we took C = 0. 
The distortions all used a triangular truncation at wavenumber 10 (T10), and computations 
were restricted to the Northern Hemisphere. Potential vorticity forecast fields were found 
to be excessively noisy near the poles, so we applied a smoother to both the forecast and 
verification fields north of 70°N. 
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Clearly the minimization greatly reduces the size of the residual error, more so for longer 
forecast periods, and most for the h500 field, less so for the sip field, and least for q315. The 
reason for this is that the distortions are limited to large scales, and the forecast errors at 
short forecast lengths, and for the q315 and sip variables, have more energy at smaller scales. 
For a qualitative assessment of the results, we turn to the forecasts of the Superstorm 93. 


4.1 Superstorm 93 

The 60 h forecast valid at 12 UTC 13 March 1993, shown alongside the verifying fields in 
Fig. 2, is at the early stage of the Superstorm 93 development. Both phase and amplitude 
errors are clearly apparent in the sea level pressure field: the center of the low is too far to 
the west nor th west, and about 30 hPa less deep than in the verifying analysis. At 500 hPa , 
the height field shows similar errors: the trough is not as deep and not as far east as in 
the analysis. The 315 K IPV field shows a southward extrusion of high IPV, which is too 
broad and weak in the forecast. The distortions applied to these forecasts (Fig. 3) all correct 
the errors to some degree, and they are generally consistent with each other, even though 
there are differences in the details. In particular, the sip corrections of the cyclone have 
a large bias correction, whereas the height and especially the IPV distortions rely more on 
displacements. As Fig. 4 shows, the original forecast errors contains more energy at the small 
scales for the IPV field, but the large-scale distortions (using the T10 truncation) manage 
to remove a large part of these errors. 

At the mature stage of development (84 hours into the forecast), the forecast captures 
the intensity of the storm, but the position of the surface low is too far to the southwest 
(Fig. 5). At 500 hPa the verification has an intense trough to the northeast of a secondary 
trough, whereas the forecast has a cutoff low in the position of the secondary trough. The 
IPV extrusion is not as wrapped around and as far to the northeast in the forecast as in 
the analysis. The distortions in this case (Fig. 6) all have northeastward displacements, and 
a pattern of bias corrections that act to weaken the surface low, and strengthen the upper 
level vorticitv (and lower the heights and pressures) to the west. 

One of the motivations for considering the IPV field was the fact that it is a quasi- 
conserved quantity. We investigated whether the distortions applied to this field would be 
less ambiguous, and exhibit more time continuity, than those applied to the height and 
pressure fields. For the cases examined so far, this holds true. An example of this is shown 
in the time series of IPV distortions for the Superstorm 93 case (Fig. 7). 


5 Covariances for background penalty terms 

We have developed an approach to determine statistics of the distortion from available data. 




SS93.060.ctrl.q315.0fg : Original forecast q315 SS93.060.ctrl.q315.0fg : Verification 







Distortion 

SS93.060.ctrl. sip. Ofg : Alignment 



Figure 3: Distortions and distorted 
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Original Forecast Error Residual Error 

SS93.060.ctrl.slp.0fg : Original forecast error sip SS93.060.ctrl.slp.0fg : Residual error 



SS93.060.ctrl.h500.0fg : Original forecast error h500 SS93.060.ctrl.h500.0fg : Residual error 


SS93.060.ctrl.q315.0fg : Original forecast error q315 SS93.060.ctrl.q315.0fg : Residual error 


■5 0 5 .5 0 5 

Figure 4: Original forecast and residual errors for forecasts at 60 h for the Superstorm 93 


case. 














Figure 6: Distortions and distorted forecasts at 84 h for the Superstorm 93 case. 
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48 h 72 h 

SS93.048.ctrl.q31 5.0fg : Alignment SS93.072.ctrl.q31 5.0fg : Alignment 



0,0 0.5 1.0 -2 0 - 1.5 -1 0 


Figure 7: Distortions at 48 h - 120 h for the Superstorm 93 case. 
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5.1 Methodology 

Parish and Derber (1992 [13]) and Parish et al. (1997 [14]) (in the JMSJ special issue), use 
a small (30 45) sample of cases of forecast divergences to develop all the statistics needed 
for the background and balance constraints. The forecast divergences are 48 — 24 h forecast 
differences for forecasts verifying at the same time. An empirical rescaling gives the 6 h 
forecast errors. This so-called “NMC” approach, has also been used successfully at EC MW F. 
At ECMYVF the empirical rescaling factor was found to be 0.9 by comparing TOVS radiance 
data to the 6 h forecast (Rabier et al. 1997 [15]). 

Using alignments to represent the forecast errors separates Jj,, the background penalty 
functional into J r , the residual penalty term measuring the misfit between the aligned back- 
ground and the analysis and Jd measuring the size of the alignment. Both of these are 
quadratic terms. J T is analogous to Jf, and would be formulated in terms of the same anal- 
ysis variables. 

The key assumptions about the forecast errors made in the NMC method are 

1. That the physical analysis variables chosen — vorticity, divergence, specific humidity, 
unbalanced temperature and unbalanced log of surface pressure are uncorrelated. 

2. That the individual spherical harmonics ( or the spherical harmonics for different total 
wavenumber n) are uncorrelated. 

3. That for each analysis variable, a vertical correlation function may be specified either 
globally or for each total wavenumber n. 

As a result of these assumptions, to calculate J(,, the model state is transformed to analysis 
variables, analyzed into spherical harmonics, and projected on the vertical EOFs of the 
vertical covariance matrices. The EOF coefficients are then squared and weighted by their 
inverse variances, which are the eigenvalues of the associated EOFs. 

Given a formulation of Jd the NMC approach could be used to derive the covariances 
for J r . We would use Jd to calculate an alignment for each forecast divergence. Then the 
residual differences would be analyzed by the NMC method. Our problem is to define the 
covariances for Jd- 

Here is a a possible approach, similar in flavor to the NMC approach. 

1. Convert forecasts to variables to be aligned on the appropriate iso-surfaces. For ex- 
ample, potential vorticity on constant potential temperature surfaces. Initially we 
examine 500 hPa data only. 

2. For each case, for each alignment variable at each level calculate an alignment without 
constraints, by using a form of J r based on a simple energy norm. Use no constraints, 
but use a forward stepwise approach to determine the spectral truncation. Initially use 
a T10 truncation. Increase the truncation in steps, stopping when the increase in the 
fraction of explained variance per degree of freedom added is less than some critical 
value (say 1%). The largest wavenumber retained will vary from case to case. Save 
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the alignment for subsequent analysis. If the alignment is saved in spectral space, be 
sure to include zero values for all wavenumbers greater than the largest wave number 
retained. Initially we examine only displacements. 

3. Develop a global vertical covariance matrix or vertical covariance matrices for each n, 
for each alignment variable from the sample of alignments. For quality control, these 
covariances should be based on trimmed samples. That is, for each alignment variable 
at each level, set the extreme values to zero. Either a critical fraction (say 3%) would 
be considered extreme, or the extreme values could be determined so that the trimmed 
sample is close to a normal distribution. 

4. The specification of Jd is now complete. It should be applied to the original sample as 
a consistency check. The alignments found with Jd should be similar to those found 
in the previous calculation of adjustments, except for cases when extreme values were 
trimmed. The residuals resulting from applying Jd can now be analyzed to obtain the 
specification of J r . 

5. As a second consistency check, calculate the alignments with the current esimates of 
J r and Jd- If the results are not substantially the same as before, go back to step 3. 

5.2 The stopping criterion 

Determining the stopping criterion for adding degrees of freedom to the representation is crit- 
ical. There are different ways of adding degrees of freedom to the distortion. In the results 
described so far, we have used the bias correction, and rotational and divergent displace- 
ments in the distortion representation, using the same T10 truncation for all 3 components. 
In the following, we describe experiments in which w r e use one (rotational displacements), 
two (rotational and divergent displacements), and all three components, using a range of 
truncation wavenumbers in each case. We have examined the results of these tests, using an 
analysis of variance approach (using an f-test stopping criterion), and an examination of the 
statistical properties of the residual error field (in particular, to what extent the correlation 
structure is homogeneous and isotropic). 

5.2.1 Analysis of Variance 

To determine a stopping criterion based on an analysis of variance, w 7 e consider the reduction 
in the residual error variance as a result of increasing the spectral truncation of the distortion 
representation. The /-statistic we compute is given by: 

_ (SSE(n 2 ) - SSEMyWfa) - d/K)) 

t ~ SS£(n 2 )/(V - d/(n 2 ) - 1) ’ 

where F* is the f-statistic corresponding to an increase of the spectral truncation from ni to 
n- 2 , SSE(n) is the residual error variance for truncation wavenumber n , df (n) is the number 
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of degrees of freedom in the corresponding distortion representation. As an estimate of the 
degrees of freedom associated with the total error field (N) we use the degrees of freedom 
contained in the spectral representation of the error field, for a rhomboidal truncation suffi- 
cient to reproduce the gridded field on the Gaussian grid used during the minimization (for 
the hemispheric fields used in these experiments, N = 1277). 

For a given truncation wavenumber ri \ , increasing the distortion truncation to n 2 results 
in a significant reduction in residual error variance if 

F* > F(p ; df(n 2 ) - d/K), N - df(n 2 ) - 1) , 


for the significance level p. 

For the potential vorticity field, we have computed the significance level of the above 
f-test for three separate cases: “default” (using all three distortion components), noamp 
(using only displacements), and “nodiv” (using only nondivergent displacements). In the 
latter case, the hemispheric mean IPV is conserved during the distortion, corresponding to a 
purely adiabatic redistribution of the original PV. For the height and pressure field, we only 
considered the “default” and “noamp” cases. In each case, the truncation of the distortion 
is a triangular truncation, and the f-test is applied to n 2 = rii + 3, i.e. for the case when the 
truncation wavenumber is increased by 3. The degrees of freedom contained in the distortion 
for the three cases is shown in Fig. 8. 



Figure 8: Degrees of freedom in the distortion, as a function of triangular truncation 

wavenumber (from T5 to T30, in increments of 3), for the hemispheric distortion for default 
(solid line), “noamp” (dotted), and “nodiv” (dashed). 

A summary of the results for the hemispheric 12-hour and 96-hour forecasts in the Su- 
perstorm 93 case is shown in Fig. 9. Results are shown for the q315, h500, and sip fields, 
for truncation wavenumbers between T5 and T25 (T30 for the “nodiv” case). The rnono- 
tonically increasing lines show the fraction of the explained variance as a function of the 
degrees of freedom in the distortion representation. For a given truncation, they are always 
highest for h500, and lowest for q315, as discussed above. It is interesting to note that, at 
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12 hours, the “default” representation results in the lowest residual error compared to the 
other representations with the same degrees of freedom, for all three fields. However, at 96 
hours, the “noamp” distortion (using displacements only) is more efficient for distorting the 
q315 field. The significance level p is unity or close to it for low truncations, but generally 
drops below 0.95 at some truncation below 7 the maximum considered in these experiments. 
The critical truncation wavenumber depends on the type of distortion, the forecast field, and 
the forecast hour. 

We have repeated these computations for the other forecast hours, and for the other 
CONTROL forecast (initialized at 00 UTC 6 March 1993), with generally similar results. 


5.2.2 Residual error correlation structure 

One of the justifications for using FCA is that the residual forecast errors are likely to be more 
homogeneous and isotropic than the original forecast errors. Visual comparison the original 
and residual error fields (Fig. 10) appear to support this hypothesis. However, the main 
difference between the two fields is the scale of variability: the residual errors have a much 
smaller scale than the original errors. This is forced by the FCA method, w 7 hich projects 
the large scale errors onto the adjustment field. To test the hypothesis of homogeneity 
and isotropy of the residual errors, we computed and compared the error correlations of 
the original and residual errors, for several forecasts, at several lead times and w'ith several 
truncations of the adjustment field. 

As before, we considered the 500 hPa height, the surface pressure and the potential 
vorticity on the 315 K potential temperature surface. The adjustment fields w 7 ere computed 
by minimizing the difference between the forecast and the corresponding analysis. It w r as 
done independently for the three variables, so that the adjustment field for the 500 hPa 
height is not the same as that for the surface pressure or the potential vorticity. 

Since this work is part of the study to determine the optimum truncation in defining the 
adjustment field, we used truncations from T05 (triangular truncation w r ith 5 waves around 
the equator) to T25. 

We first examine the 500 hPa height fields. We have the data on a regular latitude- 
longitude grid. We divide the zero to sixty degrees latitude zone into twelve equal size 
regions, sixty degrees of longitude by thirty degrees of latitude. We ignore the polar region to 
avoid having to deal w 7 ith the convergence of the meridians. Within each region we compute 
the error correlation in eight different directions and at separations of 2 to 10 degrees of great 
circle distance. To examine the homogeneity of the error fields, we average the correlations 
over all the angles within each area. Figures 11 and 12 show 7 these correlations for forecast 
starting on tw 7 o different days: 6 March 1993 and 11 March 1993. 

It can be seen from these two figures that the residual forecast error correlations behave 
in a systematic fashion as the truncation of the distortion field increases. At very low 7 
truncation the adjustment process cannot account for much of the forecast error and the 
residual errors are correlated at long distances, in a way not too dissimilar from the forecast 
error correlations. This correlation scale decreases as the truncation of the displacement 
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12 h 


ftest.SS93 12 Ctrl test: def, noamp, nodiv 
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96 h 


ftest.SS93 96 Ctrl test: def, noamp, nodiv 



o 200 400 600 800 1000 

df 
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Figure 9: Summary plot of the distortion ftest for 12 h (top) and 96 h (bottom) Superstorm 
93 forecasts. Results are shown for the q315 (black), sip (red), and h500 (blue) fields, as a 
function of the degrees of freedom retained in the distortion representation. Solid lines are 
for the “default” case, dotted for “noamp” , and dashed for “nodiv” . Shown are the fraction 
of explained variance, and the p-values of the f-test. 
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Figure 11: Error correlations averaged over all angles, in each of the 12 areas. Each plot 
shows the original error correlations in black (tropical areas) and blue (extratropical areas) 
and the residual error correlations in magenta (tropics) and green (extratropics) when the 
adjustment are computed with the indicated truncations. Twelve hour forecast from 6 March 
1993. 







- 0.2 



Figure 12: Same as Fig. 11 but for forecast from 11 March 1993 
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field increases, but little difference can be seen between the T20 and T25 truncations. At 
the same time the curves for the different regions become much more similar. This indicates 
that, beyond truncation T20 or so, the residual error correlations are much more homgeneous 
than the original forecast error correlations. It can be noticed that a large part of the 
inhomogeneity in the original error correlations is due to a large difference between the 
tropical and extra-tropical regions. This difference also appears in the residual errors, but 
to a much lesser extent. 

To evaluate the isotropy of the error correlations, we computed the correlations in eight 
different directions, every 22.5°. Figure 13 shows the spread of these curves for the original 
forecast errors and the T25 residual errors, for the 6 March 1993 forecast. The mean curves 
are averaged over all areas and all directions. The standard deviation with respect to the 
different directions is computed separately in each area, then averaged and added to or 
subtracted from the mean curve to obtain Fig. 13. It is clear that the residual error is much 
more isotropic than the original error. The figure for 11 March 1993 (not shown) is almost 
identical. 


6 Plans for future work 

A great deal of novel work has been accomplished under the current contract. In broad 
terms we have developed and tested an efficient algorithm for determining distortions. The 
algorithm and constraints are now ready for application to larger data sets, to be used to 
determine the statistics of the distortion as outlined above, and to be applied in data analysis 
by using GOES water vapor imagery to correct short term forecast errors. This future work 
is described in more detail below. A proposal (AER P778) to continue this work has been 
approved by NASA HQ in November of 1997 (UPN 622-242621), but is still under negotiation 
between AER and GSFC. An announcement of GSFC’s intent to issue a Request for Offer 
(RFO) 5-60741-253 to AER was made in CBD on 7 August 1998. 

6.1 Use of the 6.7 pm water vapor imagery 

The second set of our proposed experiments compare satellite data — in this case geosta- 
tionary 6.7 pm water vapor imagery — to a background field calculated from a short term 
forecast. The 6.7 pm water vapor imagery data are ideal for our study since they have strik- 
ing patterns and features, which can be matched by corresponding patterns and features 
in the short term forecast. Additionally geostationary water vapor data are available with 
high temporal frequency and near global coverage. However, we will begin our investigation 
with GOES data only, since the METEOSAT sensors do not have on-board calibration. (See 
Schmetz and Turpeinen (1988 [19]) for a discussion of the calibration of these data.) 

The approach for this task will be similar to that taken for the 500 hPa height fields. 
In this case, we take the short term forecast of the 6.7 pm water vapor imagery as Xj and 
the observed imagery as X a . There are two complications: the calculation of the simulated 
6.7 pm water vapor imagery, which is discussed in the next paragraph, and the need to quality 
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control the observed imagery. Quality control is required because of limited coverage, missing 
data, and the difficulty of simulating the 6.7 pm water vapor imagery at large incidence 
angles, over high terrain for dry conditions, and in the presence of cloud. It will be necessary, 
at least initially, to resample the imagery to a relatively coarse resolution. Then we will 
determine smooth displacement and amplification fields needed to best match the forecast 
and imagery. The algorithm needed here is identical to that used for the 500 hPa height 
fields. The resulting field of displacement provide a correction to the short term forecast. 

We will simulate the 6.7 pm water vapor imagery from the forecast values of temperature 
and humidity using a standard radiative tranfer model (RTM). The simulated 6.7 pm water 
vapor imagery will then be held fixed in determining the distortion. Changes in incidence 
angle related to displacements on scales of 100 km are 0(1°). The sensitivity of the calculated 
brightness temperature to incidence angle is small (Fig. 14) and will be assumed negligible 
in these calculations. The RTM used in Fig. 14 is MODTRAN (Berk et al. 1989 [2]). 
(MODTRAN is not efficient; our candidate for future calculations is OPTRAN (McMillin 
et al. 1995 [11])- During the past year we have obtained and tested OPTRAN, using a 
standard set of atmospheric profiles.) 


Simulated GOES-7 6.7 micron Brightness Temperature 



Figure 14: Variation of 6.7 pm brightness temperature with incidence angle for the U.S. 
standard atmosphere. 
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6.2 Statistics of the distortion 

To test the methodology presented in section 5 a large sample is required. Initially we wish 
to study the 500 hPa height field only. A convenient collection of forecasts and analyses for 
this purpose are the Lorenz data sets collected at ECMWF (Lorenz 1982, Simmons et al. 
1995 [10, 21]). Under the current contract we have acquired these data for four seasons— 
winter 1981, 1990, 1998; summer 1990. We plan to begin using them to develop statistics of 
the distortion as soon as the follow-on contract is finalized. 


6.3 EOF analysis 

To further limit the degrees of freedom used to represent the forecast errors, the displacement 
and amplification patterns themselves, separately and together, will be analyzed in terms of 
EOFs or rotated EOFs, to extract the typical modes of the forecast errors. The advantage of 
rotated EOFs (Richman 1981, 1986 [16, 17]), is that the resulting patterns can be localized. 
This may help to identify potential causes of model error. The EOF representation may also 
provide the basis for an enhanced version of </<*. 

We will examine the time series of the EOF coefficients of the forecast error displacement 
and amplification patterns for a fixed forecast length for correlations and periodic behavior. 
Further we will examine the evolution of forecast errors in these terms for the 1 to 5 day 
range. 


6.4 500 hPa geopotential heights assessment application 

The study of the forecast errors of the 500 hPa geopotential height fields will be extended 
to examine how forecast error varies interannually. With several years worth of data, we 
will attempt to correlate the forecast failure modes with the large-scale flow pattern and 
with other factors which might be improperly parameterized or ignored by the model. For 
example, it would be possible to represent the large scale atmospheric flow pattern and ocean 
circulation; anomalies in the SST, sea ice, and snow cover fields; and stratospheric volcanic 
aerosol by simple indices. Then these indices might be correlated with the coefficients of the 
leading EOFs of the distortion representation of the forecast errors. In this work we might 
make use of the archived forecasts and analyses of an operational model for several winters. 
However variations in predictability would then be mixed with variations in model forecast 
skill. To remove the effects of model changes it would be optimal to use a set of forecasts 
from one of the reanalysis projects. Currently reanalysis projects are underway at ECMWF, 
NMC and NASA GSFC. However, a preliminary analysis based on the Lorenz data sets is a 
logical first step. 
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